logo

This script fits state-space models to TRPD grassland bird survey data to assess species-level population trends.

Setup

Load packages & data

knitr::opts_chunk$set(echo = TRUE)
library(here)
library(png)
library(jagsUI)
library(ggplot2)
library(ggthemes)
library(tidyverse)
library(kableExtra)

source(here("Scripts", "1_Data_Import.R"))
source(here("Scripts", "2_Spatial_Data_Import.R"))
source(here("Scripts", "3_Data_Cleaning.R"))
source(here("Scripts", "4_Covariate_Data_Prep.R"))
source(here("Scripts", "5_Prep_Species.R"))

Set parameters

# Set random seed
set.seed(42563)

# Identify species and transsects for analysis.
species_to_run <- prairie_sp
transects_to_sum <- c("CHOA01", "CHOA02", "CHOA09", "CHOA14")

# ID tag for this run and create associated output directory
tag <- "20231128_plus1_to_counts"
mainDir <- here("Results", "991_bird_trends_species")
dir.create(file.path(mainDir, tag), showWarnings = TRUE)
outputDir <- here("Results", "991_bird_trends_species", tag)
subdirectories <- c("annual_n_estimates", "annual_n_plots", "diagnostic_plots", "species_level_estimates", "model_outputs", "synthesized_results", "Rhats", "model_parameters")
lapply(file.path(mainDir, tag, subdirectories), function(x) if (!dir.exists(x)) dir.create(x))

adjust_counts_by_1 <- TRUE

write.table(transects_to_sum, file = here("Results", "991_bird_trends_species", tag, "model_parameters", "transects_summed.txt"))

Refine species to run

transect_birds <- bird.tblBirds_Transects

# Merge BB and OA surveys that are actually the same transect
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB02")] <- "CHOA02"
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB09")] <- "CHOA09"
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB14")] <- "CHOA14"
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB16")] <- "CHOA14" # shouldn't be any...
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB01")] <- "CHOA01"


# Filter bird data to only specified transects
transect_birds <- transect_birds %>%
  left_join(bird.tblSurveyDefn, by = "SurveyCode") %>%
  filter(SurveyCode %in% transects_to_sum) %>%
  mutate(
    survey_string = paste0(SurveyCode, "-", SurveyDate), # unique code for survey based on area, site #, and date
    year = sub("-.*", "", SurveyDate), # convert survey date to just year
    year = as.numeric(year)
  ) # year to numeric

# Identify days and years in which all transects in the list were actually surveyed
survey_dates <- transect_birds %>%
  group_by(SurveyDate, SurveyCode) %>%
  summarize(survey_count = length(unique(SurveyDate))) %>%
  pivot_wider(names_from = SurveyCode, values_from = survey_count) %>%
  ungroup() %>%
  filter(complete.cases(.)) %>%
  select(SurveyDate)

survey_complete_years <- transect_birds %>%
  group_by(year, SurveyCode) %>%
  summarize(survey_count = length(unique(year))) %>%
  pivot_wider(names_from = SurveyCode, values_from = survey_count) %>%
  ungroup() %>%
  filter(complete.cases(.)) %>%
  select(year)

# Re-filter bird data to only dates with complete surveys across all transects. This version forces all transects to have been run on the same day.
transect_birds <- transect_birds %>%
  filter(SurveyDate %in% survey_dates$SurveyDate)

# In years surveyed more than once, pick just the last survey from that year.
surveys_all <- transect_birds %>%
  group_by(year, SurveyCode, survey_string) %>%
  summarize() %>%
  mutate(surveycode_year = paste0(SurveyCode, "-", year))
surveys_selected <- surveys_all[!duplicated(surveys_all$surveycode_year, fromLast = TRUE), ] %>%
  ungroup() %>%
  select(survey_string)

# Re-filter bird data based on previous step to only include one survey day each year.
transect_birds <- transect_birds %>%
  filter(survey_string %in% surveys_selected$survey_string)

# Only include counts > 0 (remove rows where bird only detected outside survey area).
transect_birds <- transect_birds %>%
  filter(BirdCountIn > 0)

# Get species list from data and identify species observed at least 5 times.
sp_detected <- unique(transect_birds$BirdCode)
sp_detections <- transect_birds %>%
  group_by(BirdCode) %>%
  filter(BirdCode %in% prairie_sp) %>%
  summarise(
    count_indiv = sum(BirdCountIn),
    n_surveys = n(),
    n_yrs = length(unique(year))
  )
sp_detected_above_threshold <- sp_detections %>%
  filter(count_indiv >= 5) %>%
  pull(BirdCode)

# Re-filter species to run so that it only includes species actually detected in selected surveys
species_to_run <- intersect(species_to_run, sp_detected)

# Also create version with only species detected more than 5 times.
species_to_run_above_threshold <- intersect(species_to_run, sp_detected_above_threshold)

# Save these species list for use in subsequent analyses
saveRDS(sp_detected_above_threshold, here("Results", "991_bird_trends_species", tag, "sp_detected_above_threshold.RDS"))
saveRDS(species_to_run_above_threshold, here("Results", "991_bird_trends_species", tag, "species_to_run_above_threshold.RDS"))

Format data & fit model

Define function for data prep & analysis

This code creates a function for preparing time series data for species of interest in transects of interest and then fits a state-space model (assuming density independence) to the data.

# Function to prepare time series data for species and transects of interest
prep_and_fit_CHBB <- function(species_of_interest = species_to_run, transects_to_assess = transects_to_sum) {
  transect_birds <- bird.tblBirds_Transects

  # Merge BB and OA surveys that are actually the same transect
  transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB02")] <- "CHOA02"
  transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB09")] <- "CHOA09"
  transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB14")] <- "CHOA14"
  transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB16")] <- "CHOA14" # shouldn't be any...
  transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB01")] <- "CHOA01"

  # Filter bird data to only those transects
  transect_birds <- transect_birds %>%
    left_join(bird.tblSurveyDefn, by = "SurveyCode") %>%
    filter(SurveyCode %in% transects_to_assess) %>%
    mutate(
      survey_string = paste0(SurveyCode, "-", SurveyDate), # Unique code for survey based on area, site #, and date
      year = sub("-.*", "", SurveyDate), # Convert survey date to just year
      year = as.numeric(year)
    )

  # Identify days and years in which all transects in the list were actually surveyed
  survey_dates <- transect_birds %>%
    group_by(SurveyDate, SurveyCode) %>%
    summarize(survey_count = length(unique(SurveyDate))) %>%
    pivot_wider(names_from = SurveyCode, values_from = survey_count) %>%
    ungroup() %>%
    filter(complete.cases(.)) %>%
    select(SurveyDate)

  survey_complete_years <- transect_birds %>%
    group_by(year, SurveyCode) %>%
    summarize(survey_count = length(unique(year))) %>%
    pivot_wider(names_from = SurveyCode, values_from = survey_count) %>%
    ungroup() %>%
    filter(complete.cases(.)) %>%
    select(year)

  # Re-filter bird data to only dates with complete surveys across all transects.This version forces all transects to have been run on the same day.
  transect_birds <- transect_birds %>%
    filter(SurveyDate %in% survey_dates$SurveyDate)

  # In years surveyed more than once, pick just the last survey from that year
  surveys_all <- transect_birds %>%
    group_by(year, SurveyCode, survey_string) %>%
    summarize() %>%
    mutate(surveycode_year = paste0(SurveyCode, "-", year))
  surveys_selected <- surveys_all[!duplicated(surveys_all$surveycode_year, fromLast = TRUE), ] %>%
    ungroup() %>%
    select(survey_string)

  # Re-filter bird data based on the previous step to only include one survey day each year
  transect_birds <- transect_birds %>%
    filter(survey_string %in% surveys_selected$survey_string)

  # Create df of years with complete data
  survey_years <- unique(transect_birds$year) %>%
    sort() %>%
    as.data.frame()
  survey_years <- rename(survey_years, year = .)

  # Count the number of species of interest across all transects
  sp_counts <- transect_birds %>%
    filter(BirdCode == species_of_interest) %>%
    group_by(year) %>%
    summarize(total_count = sum(BirdCountIn))

  # Zero fill for surveyed years
  sp_counts <- sp_counts %>%
    right_join(survey_years) %>%
    replace_na(list(total_count = 0)) %>%
    arrange(year)

  # This version of state-space model can't be fit when there are zeros in the data (populations that fit zero have no means of recovering/growing). Adding 1 to everything for now...
  sp_counts <- sp_counts %>%
    mutate(
      total_count_orig = total_count
    )

  if (adjust_counts_by_1 == TRUE) {
    sp_counts <- sp_counts %>%
      mutate(total_count = total_count + 1)
  }

  # Identify information about years for model
  yr_min <- min(sp_counts$year)
  yr_max <- max(sp_counts$year)
  year_range <- seq(yr_min, yr_max, 1) # Create a vector of survey years
  n_yrs <- length(year_range) # Counter for how many years of data
  year_range_df <- year_range %>%
    as.data.frame()
  year_range_df <- rename(year_range_df, year = `.`)

  # Add NAs for un-surveyed years
  sp_counts <- sp_counts %>%
    right_join(year_range_df) %>%
    arrange(year)

  #--------------------------

  # JAGS code for Bayesian hierarchical exponential population growth state space model
  {
    sink("Scripts/JAGS_code/991_bird_trends_species_version_for_MOU.jags")
    cat("
  model {

  # Specify prior values
    init_lnN ~ dnorm(ln1, 1)
    lnN[1] <- init_lnN
    r_0 ~ dnorm(0,1)  # mean intrinsic growth rate
    sigma_r ~ dunif(0, 20) # annual SD in true pop growth rate
    obs_sd ~ dunif(0,20) # annual SD in observation (measurement) error
    # convert SD to precision (tau=1/SD^2)
    sigma_tau <- pow(sigma_r, -2)
    obs_tau <- pow(obs_sd, -2)

  # Likelihood
    for (t in 1:(n_yrs-1)){
      lnN[t+1] ~ dnorm(lnN[t] + r_0, sigma_tau)
    }

  # Observation process
    for (t in 1:n_yrs) {
      lnCount[t] ~ dnorm(lnN[t], obs_tau)
      N_est[t] <- exp(lnN[t]) # convert log(N) back to real scale
    } # close t loop

  # Derived parameters

    # Geometric mean rates of change (%/year between first and last year)
    ch <- N_est[n_yrs] / N_est[1]
    trend <- (ch^(1/(n_yrs-1)))-1

    # Trend from lognormal slope
    slope_trend <- ( n_yrs * sum(1:n_yrs * lnN) - sum(1:n_yrs) * sum(lnN) ) / (n_yrs * sum((1:n_yrs)^2) - sum(1:n_yrs)^2)

  }    # End jags model
", fill = TRUE)
    sink()
  }

  # Bundle data, identify parameters to save, submit to jags
  TRPD_data <- list(
    n_yrs = n_yrs, lnCount = log(sp_counts$total_count),
    ln1 = log(sp_counts$total_count[1])
  )
  parameters <- c("r_0", "sigma_r", "obs_sd", "N_est", "trend", "slope_trend")
  TRPD1 <- jagsUI(TRPD_data,
    inits = NULL, parameters, "Scripts/JAGS_code/991_bird_trends_species_version_for_MOU.jags",
    n.adapt = 1000, n.chains = 3, n.thin = 10,
    n.iter = 300000, n.burnin = 50000, parallel = TRUE
  )

  #---------------------------------------------

  # Save model output. Models are too big to save everything, not including MCMC chains.
  mod_to_save <- TRPD1[which(sapply(TRPD1, class) != "mcmc.list")]
  saveRDS(mod_to_save, file = paste0(outputDir, "/model_outputs/", species_of_interest, "_model_outputs.RDS"))

  # Assess convergence by verifying all R-hat values < 1.1
  max.rhat <- function(x) {
    max(max_rhat <- c(sapply(x$Rhat, max, na.rm = TRUE)))
  }
  max.rhat.run <- max.rhat(TRPD1) # r.hat <1.1 suggests good convergence
  saveRDS(TRPD1$Rhat, file = paste0(outputDir, "/Rhats/", species_of_interest, "_rhats.RDS"))

  # Summarize r estimates
  r_est <- paste0("r: ", round(as.numeric(TRPD1$mean["r_0"]), 2), " (95% CI: ", round(as.numeric(TRPD1$q2.5["r_0"]), 2), "-", round(as.numeric(TRPD1$q97.5["r_0"]), 2), ")")

  # Summarize percent of r_0 posterior distribution above 0
  f_r0 <- round(as.numeric(TRPD1$f["r_0"]), 2)

  # Data frame with model estimates
  sp_values <- data.frame(
    sp = species_of_interest,
    r0_mean = TRPD1$mean$r_0,
    r0_sd = TRPD1$sd$r_0,
    r0_median = TRPD1$q50$r_0,
    r0_LCI_95 = TRPD1$q2.5$r_0,
    r0_UCI_95 = TRPD1$q97.5$r_0,
    r0_LCI_50 = TRPD1$q25$r_0,
    r0_UCI_50 = TRPD1$q75$r_0,
    f_r0 = round(as.numeric(TRPD1$f["r_0"]), 2),
    sigma_r_mean = TRPD1$mean$sigma_r,
    sigma_r_sd = TRPD1$sd$sigma_r,
    sigma_r_median = TRPD1$q50$sigma_r,
    sigma_r_LCI_95 = TRPD1$q2.5$sigma_r,
    sigma_r_UCI_95 = TRPD1$q97.5$sigma_r,
    sigma_r_LCI_50 = TRPD1$q25$sigma_r,
    sigma_r_UCI_50 = TRPD1$q75$sigma_r,
    f_sigma_r = round(as.numeric(TRPD1$f["sigma_r"]), 2),
    obs_sd_mean = TRPD1$mean$obs_sd,
    obs_sd_sd = TRPD1$sd$obs_sd,
    obs_sd_median = TRPD1$q50$obs_sd,
    obs_sd_LCI_95 = TRPD1$q2.5$obs_sd,
    obs_sd_UCI_95 = TRPD1$q97.5$obs_sd,
    obs_sd_LCI_50 = TRPD1$q25$obs_sd,
    obs_sd_UCI_50 = TRPD1$q75$obs_sd,
    f_obs_sd = round(as.numeric(TRPD1$f["obs_sd"]), 2),
    max.rhat.run = max.rhat(TRPD1),
    method = "SSMr"
  )

  write.csv(sp_values, file = paste0(outputDir, "/species_level_estimates/", species_of_interest, "_species_level_estimates.csv"))

  # Add estimates to simple dataframe
  sp_counts_w_ests <- sp_counts
  sp_counts_w_ests <- sp_counts_w_ests %>%
    mutate(
      N_est_mean = TRPD1$mean$N_est,
      N_est_sd = TRPD1$sd$N_est,
      N_est_median = TRPD1$q50$N_est,
      N_est_LCI_95 = TRPD1$q2.5$N_est,
      N_est_UCI_95 = TRPD1$q97.5$N_est,
      N_est_LCI_50 = TRPD1$q25$N_est,
      N_est_UCI_50 = TRPD1$q75$N_est,
      sp = species_of_interest,
      method = "SSMr"
    )
  write.csv(sp_counts_w_ests, file = paste0(outputDir, "/annual_n_estimates/", species_of_interest, "_annual_N_ests.csv"))

  # Plot posterior draws
  output_path <- here("Results", "991_bird_trends_species", tag, "diagnostic_plots", paste0(species_of_interest, "_diagnostic_plot.png"))
  png(output_path)
  par(mfrow = c(2, 1))

  plot(TRPD1$sims.list$obs_sd, TRPD1$sims.list$sigma_r,
    ylab = "Process error",
    xlab = "Observation error", cex = 0.3
  )
  plot(TRPD1$sims.list$r_0, TRPD1$sims.list$sigma_r,
    ylab = "Process error",
    xlab = "Process mean", cex = 0.3
  )

  dev.off()

  if (adjust_counts_by_1 == TRUE) {
    sp_counts_w_ests2 <- sp_counts_w_ests %>%
      mutate(
        N_est_LCI_95 = N_est_LCI_95 - 1,
        N_est_UCI_95 = N_est_UCI_95 - 1,
        N_est_mean = N_est_mean - 1,
        total_count = total_count - 1
      ) %>%
      mutate(
        N_est_LCI_95 = ifelse(N_est_LCI_95 < 0, 0, N_est_LCI_95),
        N_est_UCI_95 = ifelse(N_est_UCI_95 < 0, 0, N_est_UCI_95),
        N_est_mean = ifelse(N_est_mean < 0, 0, N_est_mean)
      )
  } else {
    sp_counts_w_ests2 <- sp_counts_w_ests
  }

  # Make time series plot
  plot <- ggplot(sp_counts_w_ests2, aes(x = year)) +
    geom_ribbon(aes(ymin = N_est_LCI_95, ymax = N_est_UCI_95), fill = "grey90", alpha = .6) +
    geom_line(aes(y = N_est_mean), color = "black", alpha = .7) +
    geom_point(aes(y = N_est_mean), color = "black", alpha = .7) +
    geom_point(aes(y = total_count), color = "maroon", shape = 1) +
    annotate(geom = "text", label = r_est, x = min(sp_counts_w_ests$year) + 5, y = max(sp_counts_w_ests$total_count), size = 3) +
    scale_x_continuous(breaks = seq(1960, 2021, by = 10)) +
    theme_clean() +
    labs(
      title = paste0(bird.tblBirdSpecies$BirdName[which(bird.tblBirdSpecies$BirdCode == species_of_interest)], " population dynamics"),
      subtitle = paste0("TRPD transects- ", paste(transects_to_assess, sep = "", collapse = ", "), "\nModeled estimate with 95% CI and raw counts \n", r_est),
      x = "Year",
      y = "Total count across transects"
    ) +
    theme(
      plot.title = element_text(size = 11, face = "bold"),
      plot.subtitle = element_text(size = 10, face = "italic"),
      legend.title = element_text(size = 10, face = "plain")
    )

  ggsave(plot = plot, filename = paste0(species_of_interest, "_plot.png"), device = "png", path = here("Results", "991_bird_trends_species", tag, "annual_n_plots"), width = 5, height = 3)
}

Apply function to species of interest

# # fit state space models for species of interest
# results <- lapply(species_to_run, prep_and_fit_CHBB)
# names(results) <- species_to_run
species_still_to_run <- list.files(
  path = here("Results", "991_bird_trends_species", tag, "annual_n_plots"), pattern = NULL, all.files = FALSE,
  full.names = FALSE, recursive = FALSE,
  ignore.case = FALSE, include.dirs = FALSE, no.. = FALSE
)
species_still_to_run <- gsub("\\_plot.png", "", species_still_to_run)
species_still_to_run <- unique(species_still_to_run)
species_still_to_run <- setdiff(species_to_run, species_still_to_run)

library(parallel)
numCores <- detectCores() # get the number of cores available

if (length(species_still_to_run) > 0) {
  lapply(species_still_to_run, prep_and_fit_CHBB)
  # mclapply(species_to_run, prep_and_fit_CHBB, mc.cores = numCores)
}

Visualize results

Load model results

species_level_estimates <- list.files(here("Results", "991_bird_trends_species", tag, "species_level_estimates"), pattern = "species_level_estimates.csv", full.names = TRUE) %>%
  lapply(read.csv) %>%
  bind_rows() %>%
  select(-X) %>%
  filter(sp %in% species_to_run_above_threshold) %>%
  mutate(
    r0_mean = round(r0_mean * 100, 1),
    r0_LCI_95 = round(r0_LCI_95 * 100, 1),
    r0_UCI_95 = round(r0_UCI_95 * 100, 1)
  )

Parsing annual variability

The models allow us to estimate the contribution of both observer error (undercounting/overcounting) and process error (real variability in growth rates) to the observed annual variability in populations.

error_parsing <- ggplot(data = species_level_estimates, aes(x = sigma_r_mean, y = obs_sd_mean, label = sp)) +
  geom_linerange(aes(xmin = sigma_r_LCI_50, xmax = sigma_r_UCI_50)) +
  geom_linerange(aes(ymin = obs_sd_LCI_50, ymax = obs_sd_UCI_50)) +
  geom_label(size = 2.4, label.padding = unit(.075, "lines")) +
  theme_clean() +
  labs(
    x = "Process error (Growth variate variability)",
    y = "Observer error"
  )

error_parsing

error_parsing_repel <- ggplot(data = species_level_estimates, aes(x = sigma_r_mean, y = obs_sd_mean, label = sp)) +
  ggrepel::geom_label_repel(size = 2.4, box.padding = 0.025, label.padding = .075, arrow = arrow(length = unit(0.01, "npc"), type = "open"), max.time = 20, max.iter = 100000) +
  theme_clean() +
  labs(
    x = "Process error (Growth variate variability)",
    y = "Observer error"
  )

error_parsing_repel

ggsave(plot = error_parsing, filename = "error_parsing_plot.png", device = "png", path = here("Results", "991_bird_trends_species", tag, "synthesized_results"), width = 7, height = 7)
ggsave(plot = error_parsing_repel, filename = "error_parsing_repel_plot.png", device = "png", path = here("Results", "991_bird_trends_species", tag, "synthesized_results"), width = 7, height = 7)

Annual population indices

files <- list.files(path = here("Results", "991_bird_trends_species", tag, "annual_n_plots"), full.names = TRUE)
files_names <- list.files(path = here("Results", "991_bird_trends_species", tag, "annual_n_plots"), full.names = FALSE)
files_names <- gsub("\\_plot.png", "", files_names)

files <- files[which(files_names %in% species_to_run_above_threshold)]
files_names <- files_names[which(files_names %in% species_to_run_above_threshold)]


for (i in 1:length(files)) {
  cat(paste0("#### ", bird.tblBirdSpecies.ebird$common_name[which(bird.tblBirdSpecies.ebird$BirdCode == files_names[i])], "\n"))
  cat(paste0("![](", files[i], ")\n\n"))
  cat(paste0("**Annual trend (% change/year), 1978-2023:** ", species_level_estimates$r0_mean[which(species_level_estimates$sp == files_names[i])], "% (95% CI: ", species_level_estimates$r0_LCI_95[which(species_level_estimates$sp == files_names[i])], "% to ", species_level_estimates$r0_UCI_95[which(species_level_estimates$sp == files_names[i])], "%)", "\n\n"))
}

American Crow

Annual trend (% change/year), 1978-2023: 0.9% (95% CI: -1.5% to 3.3%)

American Goldfinch

Annual trend (% change/year), 1978-2023: 2% (95% CI: -2.6% to 7.4%)

Brown-headed Cowbird

Annual trend (% change/year), 1978-2023: -0.6% (95% CI: -5.1% to 4.6%)

Blue-winged Teal

Annual trend (% change/year), 1978-2023: -0.2% (95% CI: -11.6% to 11.5%)

Clay-colored Sparrow

Annual trend (% change/year), 1978-2023: 1.1% (95% CI: -7.6% to 10.8%)

Chimney Swift

Annual trend (% change/year), 1978-2023: 0.5% (95% CI: -1.9% to 3.3%)

Common Grackle

Annual trend (% change/year), 1978-2023: -2.4% (95% CI: -10.6% to 5.1%)

Common Yellowthroat

Annual trend (% change/year), 1978-2023: 3.8% (95% CI: -7% to 15.5%)

Dickcissel

Annual trend (% change/year), 1978-2023: 1.5% (95% CI: -1.3% to 5.1%)

Eastern Bluebird

Annual trend (% change/year), 1978-2023: 2.4% (95% CI: -3.6% to 9.1%)

Eastern Kingbird

Annual trend (% change/year), 1978-2023: 1% (95% CI: -4.8% to 6.7%)

Eastern Meadowlark

Annual trend (% change/year), 1978-2023: 6.5% (95% CI: -0.8% to 13.9%)

Field Sparrow

Annual trend (% change/year), 1978-2023: 2.3% (95% CI: -9.2% to 13%)

Grasshopper Sparrow

Annual trend (% change/year), 1978-2023: 0.8% (95% CI: -4.7% to 6.9%)

Henslow’s Sparrow

Annual trend (% change/year), 1978-2023: 5.2% (95% CI: -6.8% to 16.9%)

Killdeer

Annual trend (% change/year), 1978-2023: 1% (95% CI: -3.6% to 5.8%)

Mourning Dove

Annual trend (% change/year), 1978-2023: -2.3% (95% CI: -10% to 5.6%)

Northern Cardinal

Annual trend (% change/year), 1978-2023: 0.5% (95% CI: -5.4% to 6.2%)

Ring-necked Pheasant

Annual trend (% change/year), 1978-2023: 0.9% (95% CI: -1.5% to 3.9%)

Red-winged Blackbird

Annual trend (% change/year), 1978-2023: 1.2% (95% CI: -11.2% to 13.9%)

Savannah Sparrow

Annual trend (% change/year), 1978-2023: -1.2% (95% CI: -7.7% to 5.1%)

Sedge Wren

Annual trend (% change/year), 1978-2023: 2.7% (95% CI: -5.5% to 11.6%)

Song Sparrow

Annual trend (% change/year), 1978-2023: -1.4% (95% CI: -13.9% to 10.9%)

Tree Swallow

Annual trend (% change/year), 1978-2023: 5.7% (95% CI: -3.9% to 15.4%)

Vesper Sparrow

Annual trend (% change/year), 1978-2023: -2.2% (95% CI: -9.1% to 6.5%)

Western Meadowlark

Annual trend (% change/year), 1978-2023: -1.7% (95% CI: -9.6% to 5.5%)

Willow Flycatcher

Annual trend (% change/year), 1978-2023: 0.3% (95% CI: -4.6% to 4.8%)


Print session info:

sessionInfo()
## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22621)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] jagsUI_1.5.2           png_0.1-8              transformr_0.1.3       gifski_1.12.0-2       
##  [5] gganimate_1.0.8        ggspatial_1.1.9        ggrepel_0.9.4          RVAideMemoire_0.9-83-7
##  [9] pairwiseAdonis_0.4.1   cluster_2.1.4          goeveg_0.6.5           vegan_2.6-4           
## [13] lattice_0.20-45        permute_0.9-7          mapview_2.11.2         ggtext_0.1.2          
## [17] ratelimitr_0.4.1       rvest_1.0.3            trelliscopejs_0.2.6    plotly_4.10.3         
## [21] auk_0.7.0              readxl_1.4.3           kableExtra_1.3.4       ggthemes_5.0.0        
## [25] forcats_1.0.0          stringr_1.5.1          purrr_1.0.2            readr_2.1.4           
## [29] tidyr_1.3.0            tibble_3.2.1           ggplot2_3.4.4          tidyverse_2.0.0       
## [33] lubridate_1.9.3        dplyr_1.1.4            DT_0.31                sf_1.0-14             
## [37] rgdal_1.6-4            sp_2.1-2               RODBC_1.3-23           here_1.0.1            
## 
## loaded via a namespace (and not attached):
##   [1] spam_2.10-0             uuid_1.1-1              backports_1.4.1        
##   [4] Hmisc_5.1-1             systemfonts_1.0.5       plyr_1.8.9             
##   [7] lazyeval_0.2.2          splines_4.2.2           crosstalk_1.2.1        
##  [10] leaflet_2.2.1           unmarked_1.3.2          leafpop_0.1.0          
##  [13] digest_0.6.33           htmltools_0.5.7         leaflet.providers_2.0.0
##  [16] fansi_1.0.5             magrittr_2.0.3          checkmate_2.3.0        
##  [19] tzdb_0.4.0              vroom_1.6.4             svglite_2.1.2          
##  [22] timechange_0.2.0        lpSolve_5.6.19          prettyunits_1.2.0      
##  [25] ggh4x_0.2.6             colorspace_2.1-0        leafem_0.2.3           
##  [28] textshaping_0.3.7       xfun_0.41               crayon_1.5.2           
##  [31] jsonlite_1.8.7          brew_1.0-8              glue_1.6.2             
##  [34] gtable_0.3.4            webshot_0.5.5           maps_3.4.1.1           
##  [37] scales_1.3.0            DBI_1.1.3               Rcpp_1.0.11            
##  [40] htmlTable_2.4.2         viridisLite_0.4.2       xtable_1.8-4           
##  [43] progress_1.2.3          gridtext_0.1.5          units_0.8-5            
##  [46] foreign_0.8-83          bit_4.0.5               proxy_0.4-27           
##  [49] mclust_6.0.1            dotCall64_1.1-1         Formula_1.2-5          
##  [52] stats4_4.2.2            htmlwidgets_1.6.4       epuRate_0.1            
##  [55] httr_1.4.7              wk_0.9.1                ellipsis_0.3.2         
##  [58] pkgconfig_2.0.3         farver_2.1.1            nnet_7.3-18            
##  [61] sass_0.4.7              utf8_1.2.4              tidyselect_1.2.0       
##  [64] labeling_0.4.3          rlang_1.1.2             later_1.3.1            
##  [67] munsell_0.5.0           cellranger_1.1.0        tools_4.2.2            
##  [70] cachem_1.0.8            cli_3.6.1               generics_0.1.3         
##  [73] evaluate_0.23           fastmap_1.1.1           yaml_2.3.7             
##  [76] ragg_1.2.6              knitr_1.45              bit64_4.0.5            
##  [79] s2_1.1.4                satellite_1.0.4         pbapply_1.7-2          
##  [82] nlme_3.1-160            mime_0.12               xml2_1.3.5             
##  [85] compiler_4.2.2          rstudioapi_0.15.0       e1071_1.7-13           
##  [88] tweenr_2.0.2            bslib_0.6.1             stringi_1.8.2          
##  [91] highr_0.10              fields_15.2             Matrix_1.6-4           
##  [94] classInt_0.4-10         commonmark_1.9.0        markdown_1.12          
##  [97] vctrs_0.6.5             pillar_1.9.0            lifecycle_1.0.4        
## [100] jquerylib_0.1.4         MCMCvis_0.16.3          data.table_1.14.8      
## [103] cowplot_1.1.1           raster_3.6-26           httpuv_1.6.12          
## [106] R6_2.5.1                promises_1.2.1          gridExtra_2.3          
## [109] KernSmooth_2.23-20      rjags_4-15              codetools_0.2-18       
## [112] MASS_7.3-58.1           assertthat_0.2.1        rprojroot_2.0.4        
## [115] withr_2.5.2             autocogs_0.1.4          mgcv_1.8-41            
## [118] hms_1.1.3               terra_1.7-55            rpart_4.1.19           
## [121] grid_4.2.2              rebird_1.3.0            coda_0.19-4            
## [124] class_7.3-20            rmarkdown_2.25          DistributionUtils_0.6-1
## [127] shiny_1.8.0             base64enc_0.1-3
 




By Sam Safran